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Abstract 

O 

■ We investigate the crystal structure of classical systems of spherical par- 
ticles with an embedded point dipole at T = 0. The ferroelectric ground state 
energy is calculated using generalizations of the Ewald summation technique. 
Due to the reduced symmetry compared to the nonpolar case the crystals 
are never strictly cubic. For the Stockmayer (i.e., Lennard- Jones plus dipo- 
lar) interaction three phases are found upon increasing the dipole moment: 
hexagonal, body-centered orthorhombic, and body-centered tetragonal. An 
even richer phase diagram arises for dipolar soft spheres with a purely re- 
pulsive inverse power law potential ~ r~ n . A crossover between qualitatively 

£T) | different sequences of phases occurs near the exponent n = 12. The results are 

applicable to electro- and magnetorheological fluids. In addition to the exact 
q ■ ground state analysis we study freezing of the Stockmayer fluid by density- 

■ functional theory. 

PACS numbers: 61.50.Ah, 83.80.Gv, 77.80.-e, 64.70.Kb 
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I. INTRODUCTION 



The last years have seen a revival of interest in simple dipolar fluids, which consist of 
spherical particles with embedded point dipoles, triggered by two unexpected observations: 
At high interaction strengths and high densities a liquid phase with long-range ferroelectric 
orientational order but without any positional order occurs |J; at low densities the particles 
form a gas of chains behaving like living polymers []2| . The latter effect was suggested as an 
explanation for the apparent absence of gas-liquid condensation in dipolar hard spheres [[3] , 
but this subject is still under discussion ||. Both phenomena have first been detected in 
computer simulations followed by theoretical work [[TT]|-|H|] . This discussion applies to 

electric and magnetic dipoles in complete analogy; in the following we will use the electric 
language. 

Knowledge about the solid phase of these systems is necessary in order to know under 
which circumstances formation of the ferroelectric liquid is preempted by freezing. This ques- 
tion has been tackled theoretically within two different versions of density-functional theory. 
Groh and Dietrich Jl5| found a stable ferroelectric liquid for the Stockmayer (i.e., dipolar 
Lennard- Jones) model if the dipole moment is high enough. Klapp and coworkers, however, 
found within their approach that the ferroelectric liquid is always metastable in comparison 
with the solid both for Stockmayer and dipolar hard sphere fluids [^,|T^]. But in a recent 
study they demonstrated that their result depends sensitively on the applied approxima- 
tions |TJ|. The only solid structures considered in these papers |TB|-|TBJ are face-centered 
cubic (fee), which is the known crystal structure in the non-polar limit, and body-centered 
tetragonal (bet) with the special axis ratio c/a = \/2/3, which had been determined before 
as the ground state of dipolar hard spheres [IS] and has also been observed in simulations 
20| and in experiments |2T|j22]] with electrorheological fluids, i.e., suspensions of polarizable 



colloidal particles. 

In a recent simulation Gao and Zeng |23] confirmed the occurrence of a stable ferroelectric 
liquid phase for the Stockmayer model and, in addition, determined for the first time portions 
of the phase boundaries between isotropic and ferroelectric liquid and between ferroelectric 
liquid and solid. Although these results support the basic conclusion in Ref. JDJ they found 
that the stable solid phase of the Stockmayer model is body-centered orthorhombic (bco) 
and additionally observed a metastable distorted hexagonal structure, possibilities that have 
not been taken into account in the theoretical work f!5|-|18| so far. Both bco and bet as well 
as an fee crytal with helically varying polarization direction have been reported before in 
simulations of dipolar hard spheres pi] , but the thermodynamically stable state could not 
be determined. Certain structures may be suppressed in simulations due to the periodic 
boundary conditions if the cell shape is not flexible enough. 

A simple heurestic argument why a cubic structure is not expected runs as follows. All 
these crystals are ferroelectric and hence have less symmetry than in the nonpolar case; 
the point symmetries can only be reflections at planes that contain the polarization axis 
and rotations around this axis. Therefore if, e.g., in a cubic crystal a polarization in [100] 
direction is switched on, the additional introduction of a contraction along this direction 
does not further reduce the remaining symmetry. Hence generically the crystal will have an 
axis ratio different from unity, i.e., it will be tetragonal. For polarization along the initial 
[110] and [111] directions the same reasoning leads to orthorhombic and trigonal crystals. 
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Thus a ferroelectric solid can never be strictly cubic, in contrast to the assumptions made 
in the density-functional work described above and in similar work on the Heisenberg fluid 

Thus up to now it is quite unclear which crystal structure(s) actually are stable in 
these widely used dipolar model systems. As a first step to elucidate this region of their 
phase diagram in the present work we determine the ground state structure as a function of 
density, dipole strength, and softness of the isotropic interaction potential. Since in principle 
an infinite number of crystal structures exists, characterized by an increasing number of 
parameters with increasing number of particles in the basis, an exhaustive search for the 
ground state is not possible. However, we include a rather large class of candidate structures, 
comprising among others all reasonable simple Bravais lattices. A relatively complex phase 
behavior is found at T = 0. There is reason to expect that the ground state results remain 
qualitatively valid at not too high temperatures and thus provide a valuable starting point 
for a more complete analysis of the full phase diagram. 

In a recent work |26[ ground state energies and ordering temperatures for ferroelectric 
and antiferroelectric arrangements of Ising dipoles on several lattices have been determined, 
but no attempt has been made to optimize the lattice structure. In this study also the 
presence of isotropic interactions was not taken into account. 

In our analysis we always assume a spatially homogeneous polarization throughout the 
sample. In real ferroelectric materials domains will form due to the long range nature of the 
dipolar interaction [27] . The domain structure in the liquid ferroelectric phase for a cubic 
sample shape has been analyzed in detail in Ref. p8| . This more complicated situation is 
avoided in two cases: (i) a needle-shaped sample (infinite aspect ratio), implying a vanishing 
depolarization factor; (ii) cancelling of the induced surface charges by free charges which are 
present in a conducting surrounding medium or as impurities in the material itself. 

After having reached a basically complete overview of the thermodynamically stable 
ground states in such systems, which is interesting in its own right, in a second step by using 
density-functional theory we return to the issue which solid phase the ferroelectric liquid 
phase the Stockmayer fluid forms upon freezing. This closes the aforementioned gap between 
the analysis of the solid phases as considered so far in the theoretical analyses |I5}{18[] and 
the observation of those types of solid structures as found in simulations |23| . It consolidates 
the theoretical prediction |15j that the Stockmayer fluid can exhibit a thermodynamically 
stable ferroelectric liquid phase, in accordance with the simulation results. 



II. MODELS AND METHODS 

We study systems of spherical dipolar particles interacting via the dipolar potential 

. . — m 2 r 2 + 3(m • r) 2 , , 

w dip (r) = (1) 

where r is the interparticle vector and m the dipole moment, which is assumed to have the 
same orientation for all particles. In addition there is an isotropic interaction Wi SO which is 
taken as either the purely repulsive soft sphere (SS) potential 

w SS (r) = 4e (-J (2) 
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or the Lennard- Jones (LJ) potential 



w L j{r) = 4e 



r 



(3) 



The parameters e and a define the energy and length scales, respectively. The ground state 
properties only depend on the reduced dipole moment m* = (m 2 /(cr 3 e)) 1//2 and the reduced 
particle density p* = a 3 N/V of iV particles in the volume V. Since in the SS case only 
the combination ea n enters the potential, the only independent thermodynamic parameter 
is x — m 

2/ p *n/3-l_ J n the L J 

case the ground state energy per particle U has the form 
U = eE(p*,m*) whereas for the SS case one has U = ep* n ^E(m* 2 / p* n / 3_1 ) where E and 
E are dimensionless scaling functions which contain all the structure dependence. For soft 
spheres the optimum structure remains the same as long as x is kept constant because 
the prefactor 

p *n/3 ig 

the same for all structures belonging to a given reduced density p*. 
Furthermore one can show that also phase coexistence densities p* determined for one value 
of m* can be scaled to another dipole moment m*' according to p*\ = (m*/m*') 6//< - n_3 )p*, 
because d(p*U)/dp* = ep* n ^ 3 Ei(x) with another scaling function E\. Hence for the SS model 
without loss of generality we set m* — 1. 

Explicitly the ground state energy per particle is 

U = U dip + U iso = [w dip (R + T) +w iso (R + r)] (4) 

R,T 

where R runs over the Ni lattice vectors of a Bravais lattice and r over the M positions 
of the basis particles within one unit cell so that N = iVjM; the prime on the summation 
sign indicates that the term with R + t = must be omitted. Here we implicitly have 
replaced a double sum over the lattice sites by a single sum, assuming that the average 
energy per particle is equal to the energy of the particle at the origin. It is a non-trivial issue 
whether this assumption is justified for the long-ranged dipolar potential. In the appendix 
we show that it is correct for a homogeneous orientational configuration in an ellipsoidal 
sample shape, but not, e.g., for a parallelepiped. (But we recall that only under the conditions 
mentioned in the last but one paragraph of Sec. | the ground state has spatially homogeneous 
orientational order as assumed in Eq. (|4]).) Straightforward numerical calculation of these 
lattice sums is hampered by their slow convergence. Therefore generalizations of the Ewald 
technique are employed to transform the sums into a more rapidly converging form. The 
basic idea for evaluating adroitly a general sum 5^ R /(R) is to write /(R) = h(R) + g(R) 
where h(R) decays rapidly in real space and the Fourier transform ^(k) = J d 3 r e _lk ' r g(r) 



decays rapidly in reciprocal space ||29|| . The sum of g is then evaluated in reciprocal space 
using the Poisson sum formula. For the inverse power sums with /(R) = R~ n one uses 
h(R) = r(n/2,i> 2 R 2 )/(R n r(n/2)) where T(a,x) is the incomplete Gamma function and 
obtains for a Bravais lattice 



R^O v ' ' R^O k^O v 7 



2 „ , 2 vr 3 / 2 n _ 3 

— v H v 

n n — 3 V c 



(5) 
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Here k runs over the reciprocal lattice, V c is the volume of the unit cell and the last two 
contributions take into account the omissions of the terms R = and k = 0. The parameter 
v can be chosen arbitrarily and the independence of v of the total sum provides a convenient 
check of the algorithm. In practice v is chosen such that both sums converge approximately 
with the same rate. Typically a few hundred lattice vectors in real and reciprocal space 
are sufficient to obtain machine precision (10 -16 ). The straightforward generalization to M 
particles in the basis leads to 



U. 
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• (6) 



Clearly for the L J case C/j SO = £/#5 — t/55 • 



An analogous expression for the dipolar lattice sum can be obtained from the well-known 



result for the potential of a lattice of point charges [31 
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with 



^(R) = ^l(R 2 - 3^ - 2Ml®e-'>* + ^erfc^). 



(8) 



Due to the slow decay of the dipolar sum is actually only conditionally convergent, 
i.e., the result is shape dependent. In the appendix we present the derivation of Eq. (|7|) and 
show that it corresponds to the case of a needle-shaped sample which is of interest here. The 
corresponding generalization to lattices with a basis reads 



dip 
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(9) 



If the dipoles are not oriented along the z axis of the lattice but have a general orientation 
m the dipole sum has the form 'Y^ij^iTijifij with a symmetrical matrix (see Eq. ([!])). 
The optimum direction is then necessarily along one of the eigenvectors of T which coincide 
with high symmetry lattice directions. Thus it suffices to consider only one or two possible 
orientations for each lattice. 

The following crystal structures, displayed in Figs. [I] and ^|, were included in the search 
for the minimum of the energy (except for trigonal lattices the polarization is always along 
the c axis): 

1. body-centered orthorhombic (bco) with axis lengths a, b, c; reduces to bet for a = b, 
to fee polarized along [110] for b/a = c/a = l/y/2 (see Fig. [I]), and to fee polarized 
along [001] for b/a = 1, c/a = y/2. 
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2. 



face-centered orthorhombic (fco); note that for the tetragonal case (6 = a) face-centered 
and body-centered lattices are equivalent. 



3. trigonal (trig): three equal axes with angle 7 between any pair of them, polarized along 
[111]; reduces to various cubic lattices polarized along [111] for special values of 7. 

4. hexagonal with axis lengths a, c and a second basis particle at r = a/2(l, 1/V3, c/a) 
polarized along the c axis (hexc); corresponds to hexagonal close packed (hep) for 
c/a = ^8/3. 

5. an orthorhombic lattice with four basis particles at Tq = 0, T\ = (a/6, 6/2, c/2), T2 = 
(a/2, 0,c/2), and r 3 = (2a/3, 6/2,0) which can be viewed as a distorted hep lattice 
with polarization in the a6-plane (hexab); this structure was observed in Ref. |23| . 



These possibilities have been chosen for the following reasons. The structures bco, hexc, 
and hexab are generated by slight distortions of the close packed fee and hep lattices which 
represent the ground state for m = 0. The structures fco and trig are included, because they 
approach reasonable low density configurations in certain limits, which have been observed in 



electrorheological fluids [22]. For 7 — > 27r/3 trig degenerates into a hexagonal array of dipolar 
chains, and for c,6>o fco develops into a collection of parallel sheets. The only remaining 
Bravais lattices, monoclinic and triclinic, are improbable and also difficult to handle because 
of the larger number of free parameters. We tested a monoclinic variation of bco and always 
found that the energy is minimized for a right angle between the axes. Concerning more 
complex lattices it is difficult to define a reasonable parameter space without a physically 
motivated structure to start from. 



III. RESULTS 

A. Stockmayer model at T = 

If one starts from the fee structure of the nonpolar LJ system and introduces a dipole 
moment, all orientations of the polarization have the same energy, as long as the cubic sym- 
metry is preserved. However, as discussed in the introduction, the crystal actually distorts 
towards bco, bet, or trig, depending on the direction of the polarization. The numerical find- 
ings indicate that the [110] direction is selected, corresponding to a bco distortion as shown 
in Fig. |l|. Intuitively this can be understood as a preference for the formation of chains along 
the polarization direction: [110] points towards the nearest neighbors so that the particle 
distance along the chains is smallest in this case. Indeed, upon increasing the dipole mo- 
ment the length of the c axis decreases, reflecting the strong attractive interactions along the 
chains. Figure ^| displays the dependence of the bco axis ratios on the reduced dipole moment 
m*. The result of Gao and Zeng |23J for p* = 1.24, m* = 2.5, and T* = k B T/e = 0.7 denoted 
by the diamonds lies very close to the T = result. Obviously the ground state gives a good 
approximation to the equilibrium state at least up to half the triple temperature T t * ~ 1.3 
|2~3fl . The same authors also report a metastable hexab structure at m* = 2.5, T* = 0.8, and 



p* = 1.146 with the axis ratios c/a = 0.497 and b/a = 0.953. The corresponding values in 
the ideal hep lattice are c/a = 0.577 and b/a = 0.942, respectively. At T = the hexab 
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crystal is only metastable, too, with c/a = 0.510 and b/a = 0.961. With increasing dipole 
moment the two bco axis lengths b and a perpendicular to m approach each other until at 
a critical value m* a continuous transition to bet takes place (Fig. |3]). After a cusp at m* in 
the bet phase c/a decreases again and is much lower than the "ideal" value a/2/3 = 0.816 
found for dipolar hard spheres. 

The true ground state of the nonpolar LJ model, however, is hep which has an energy 
very slightly below the fee value (by 0.01-0.02% depending on density). This difference is 
only due to second nearest neighbors because the number and distances of the 12 nearest 
neighbors are equal in fee and hep. The dipolar energy of the hep lattice with polarization 
along the c axis |20[ is lower than with polarization along the a axis |19[ . Therefore for small 
dipole moments a slightly contracted hep lattice (hexc) is preferred over bco, although the 
energy differences are always very small, e.g., below 0.16% for p* = 1.24. At larger dipole 
moments bco takes over as the stable phase. Because of the smallness of these differences 
entropy effects at finite temperature may easily tip the balance towards one or the other 
phase. The full T = phase diagram is shown in Fig. |j. The hexc-bco transition is first 
order but exhibits only a tiny density jump Ap* ~ 5 • 10~ 3 . The dashed line Pmin( m *) 
marks the absolute minimum of the energy per particle over all densities. If a system is 
prepared with a density below p* min it spontaneously shrinks to this minimum and leaves a 
corresponding portion of empty space. Thus the region p* < p* min is a two-phase coexistence 
region between the lowest energy solid and an infinitely diluted gas. For T > it connects 
to gas-solid coexistence while the liquid phase (s) appear only at higher temperature above 
a triple point. That p* min corresponds to a phase coexistence density can also be shown 
more formally by performing a double tangent construction on the (free) energy density 
u(p) = pU(p). For finite temperatures an entropic term ~ Tplnp must be added for the gas 
phase, so that a double tangent at p„ and p s can be constructed with 



0. 



du 
dp 



1 1 



P 



for T -> 0. 



(10) 



The last equation is equivalent to the minimum condition dll/dp = 0. 

Thus we conclude that at T = the Stockmayer crystal in the commonly studied range 
m* < 3 is either hexc or bco, but neither fee nor bet as assumed in various studies before. 



B. Stockmayer model at T > 



In this section we present the predictions of the density- functional theory (DFT), which 
we applied to freezing of the Stockmayer fluid in our previous work ||15|| , when the addi- 
tional crystal structures hexc and bco are taken into account. This approach is based on 
a perturbation expansion around the hard sphere solid which is treated in the modified 



weighted-density approximation of Denton and Ashcroft |32| . The long-ranged isotropic and 
dipolar interactions are added in such a way that a successful theory of the Lennard- Jones 
fluid |$3| is reproduced in the nonpolar limit. The detailed definition of the density functional 
is given in Ref. |l5j and therefore is not repeated here. 

In order to treat the hexc phase with more than one basis particle some of the expressions 
given in Ref. must be generalized; e.g., the expression for the long-range contribution 
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to the excess free energy turns into 



^p- = -|A. 2 «? ( 1 - 6£|S(k)| 2 P 2 (cos^)e^M^ ) (11) 



with S(k) = M 1 ^ T exp(— • r) and the other quantities as defined in Ref. |15| . Besides 



the peak width 7 and the orientational order parameters, the axis ratios now appear as 
additional variables in the minimization. 

Figure [5] shows the variation of the bco axis ratios as function of the dipole moment at 
fixed temperature and density. For low m* the particles in the crystal are orientationally 
disordered and both axis ratios are equal to I/a/2 corresponding to an fee lattice. When 
ferroelectric order sets in both axis ratios eventually decrease but exhibit a peculiar mini- 
mum and maximum in an intermediate range of values for m*. The contraction along the 
polarization direction (c axis) is in accordance with the ground state result. On the other 
hand, in contrast to the DFT results, both at T = and in the simulations the value of 
b/a is larger than l/y/2 and increases with m* (compare Fig. ||). A quantitative comparison 
with the simulation results for m* = 2.5 is not possible because such high dipole moments 
could not be reached in the theory. (For large values of m* the crystalline density peaks 
become very narrow which leads to convergence problems for example in Eq. flll]).) In DFT 
usually both axis ratios decrease with increasing dipole moment, decreasing temperature, or 
decreasing density. 

For m* = 2 hexc turns out to be the stable crystal structure at all temperatures. In 
Fig. |6| we show the calculated ferroelectric liquid-ferroelectric solid and gas-ferroelectric 
solid transition densities. For comparison our previous data assuming an fee structure are 
displayed, too. The shift due to the new crystal structure is relatively small and has no 
impact on the occurrence of the ferroelectric liquid phase as such. The axis ratio c/a always 
lies below the ideal hep value a/8/3 = 1.633; it varies between 1.58 and 1.62 along the parts 
of the coexistence lines shown in Fig. ||. Surprisingly it decreases with increasing density 
while the opposite trend is observed in the ground state. For the lower dipole moment 
m* = 1 the free energy differences between hexc, bco, and fee are smaller than the numerical 
accuracy so that with the present tools it is not possible to decide which phase is the stable 
one. In any case, however, the loci of the phase boundaries remain practically the same as 
calculated in Ref. Hl5| . Thus we conclude that our DFT prediction of the occurrence of a 



stable ferroelectric liquid phase within a certain parameter range is not altered by taking 
into account more possibilities for the crystal structure and thus is in agreement with the 
findings of the simulations. 



C. Dipolar soft spheres at T = 

In this case fee is more stable than hep for m* = 0, or, equivalently, for m* 7^ and 
p* — > 00. For large exponents n with decreasing density qualitatively the same sequence of 
transitions fcc-bco-bct occurs as discussed for the Stockmayer system. But as shown by the 
phase diagram in Fig. [7] within the bco range hexc is stable in an intermediate range. Also 
in contrast to the Stockmayer model the vapor phase can coexist only with the bet solid. 
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For n — > oo one recovers the limit of dipolar hard spheres. Here the lowest energy state 
is the "ideal" bet with c/a = a/2/3 at p* = 4/3 Jl9| which coexists with an infinitely diluted 
state at p* = 0. For higher densities close packed bco structures with 

b ---( C -Y (12) 
a p* \aJ 

occur in which each particle has ten nearest neighbors at distance o. At p* = 1.383 a two- 
phase bco-hexc region starts that extends up to the maximum possible density p* = \[2 
where the ideal hep lattice with c/a = a/8/3 is stable. For n — > oo the values of the axis 
ratios along the various phase boundaries in Fig. [7] converge towards the aforementioned 
hard sphere values. 

A quite different behavior occurs for small exponents n. Here the isotropic repulsion 
cannot be overcome by the dipolar attraction in a three-dimensional structure so that the 
energy is lowest at p = which means that the solid phase trig remains stable down to 
arbitrarily low densities without encoutering a gas-solid phase separation. 

A semi-quantitative understanding of this effect can be obtained by considering arrange- 
ments of parallel chains, to which all crystal structures degenerate for p — > 0. The intrachain 
energy per particle of a soft sphere chain is |19| 
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777 ^ / fj\7l 

U ch = -2((3) — + 4e(-) C(n) (13) 

where ((n) = Y^kLi k~ n is the Riemann zeta function and a the particle distance. The 
equilibrium distance follows by minimization: 

Oeq = f 2n((n) 
" V3C(3)m* 2 



a 



(14) 

The purely dipolar interaction energy between two parallel chains with distance r and lon- 



gitudinal offset z was calculated by Tao and Sun [19 : 



T jdip 167r 2 m 2 / 2nkr\ ( 

■\^ e " 2Wacos (^?) 



a 

, , r- (15) 



where K v denotes modified Bessel functions. Using similar methods one finds for the isotropic 
contribution of a single power-law repulsion 



rriso Ua " 



ch— ch 



a r(n/2) 



r ( ! T i ) ,f /MX,, /27rA;r\ /2tt^\ 



(16) 



All but the first term decay exponentially for r — > oo. The total chain-chain interaction 
f^ch-ch = ^h-ch + ^ch-ch' evaluated at a = a eq and z = a/2 behaves qualitatively different 
for large and small exponents n, as shown in Fig. II. For large n it exhibits a minimum with 
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£4h-ch < at intermediate r followed by a maximum and eventually an algebraic decay for 
r — > oo, whereas for small n it is repulsive and monotonously decreasing for all distances 
r. The minimum reaches zero at n = 11.93 and disappears at n — 11.89. Thus the most 
frequently used value 12 for the exponent, mainly chosen for historical reasons, just marks 
the crossover between attractive and repulsive soft dipolar chains which is also reflected in 
the phase diagram. Near n = 12 a low density fco phase appears. In this phase for p — > the 
particles form parallel sheets, due to the small attraction between chains (see Fig. [Sp. For 
slightly smaller values of n the ground state becomes trigonal with opening angle 7 — > 2tx/2> 
for p — > 0, i.e., a hexagonal lattice of chains with relative longitudinal shifts ±a/3 between 
nearest neighbors. Using the equations above one can verify that this limiting structure is 
more favorable than a hexagonal arrangement with shifts and a/2 which can be obtained 
from a bco lattice with b/a = a/3. In the same region the bet phase disappears so that the 
transition sequence upon increasing density becomes trig-bco-hexc-bco-fcc. 

In both models the hexab phase turns out to be metastable for all parameters. All solid 
lines are first order transitions. Except for the transitions bco-hexc near n = 00 and fco-bct 
near n = 12 the density gaps are always very small. Dashed lines denote the continuous 
bco-bct and fco-bct transitions. 



IV. DISCUSSION 

Even at T = the investigated model systems show a rich phase behavior with a variety 
of solid-solid phase transitions as function of dipole moment, density, and softness of the 
repulsion (see Figs. |] and |7]). Concerning the experimental relevance of these phenomena 
molecular dipolar fluids are natural first candidates. However, for them quantitative com- 
parison is impeded by the fact that typically such particles exhibit additional short-ranged 



steric anisotropics which become important at freezing densities |34|]. Such a sensitive depen- 
dence of the phase behavior on the details of the repulsive part of the interaction potential is 
supported by our findings concerning the decay exponent n (see Fig. |7]). The closest effective 
realization of our models are colloidal suspensions of monodisperse spherical particles. The 
dipole moment can either be a permant magnetic moment as in ferrofluids |35[ or an induced 



electric or magnetic moment as in electrorheological (ER) or magnetorheological fluids |36 
At present it is still difficult to prepare stable ferrofluids at sufficiently high densities. How- 
ever, recently the occurrence of gas-liquid-solid phase transitions in nearly monodisperse 
solutions of maghemite 7-Fe203 nanoparticles in water has been reported [J37| . In this ionic 
ferrofluid the effective isotropic interaction between the particles can be tuned by changing 
the screened Coulomb interaction via adding salt. This kind of intervention into the interac- 
tion potential is very interesting since our analysis demonstrates that, as mentioned above, 
the occurrence of different solid phases depends sensitively on the details of the isotropic 
repulsion. 

The electrostatic energy of an arrangement of polarizable ER spheres with radius a and 



dielectric constant ep in a solvent of dielectric constant ep and an external field E is |L9] 



ae F a 3 E 2 

ER ~ 2(l + 2aa»0y 1 } 
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with a = (ep — ep)/(ep + 2ep) and = Udi P /m 2 the corresponding reduced energy for 
permanent dipole moments m. For small a (\a\ < 1/2 by definition, a ~ 0.3 has been es- 
timated for a silicon oil ER fluid P5| , while a > —1/2 for water based fluids) Eq. ( [T7| ) can 
be expanded and the structure dependent terms take on the same form as for permanent 
dipoles with an effective moment = a 2 epa 6 E 2 so that the calculated phase diagrams 
apply without changes. Moreover with typical values [^] E = lkV/mm and a = 0.5/mi 
the dipolar energy at contact m 2 s /(2a) 3 is larger than the thermal energy ksT at room 
temperature by three orders of magnitude, justifying the use of our ground state analysis. 
Concerning the effective isotropic interactions, dispersion forces as modelled by the LJ po- 



tential are present in ER fluids |39[ but usually neglegibly small ||40|| . Steric repulsion at 
small distances is achieved by polymer coating. The length and density of polymers deter- 
mines the softness of the repulsion although it will be difficult to reproduce a power law 
dependence as considered theoretically above. Thus chemical tayloring of the particle surface 
represents an option to produce softly repulsive potentials which differ from the hard sphere 
behavior usually assumed in ER models. We expect that our calculated phase diagrams 
at least qualitatively reflect the behavior of such ER fluids. While the subtle issue of the 
relative stability of the fee and hep phases and the corresponding bco-hexc transitions are 
probably masked by neglected effects such as higher multipoles, many-particle interactions, 
and polydispersity, the phase sequence fcc-bco-bct upon increasing field strength should be 
insensitive to these details. 

At small field strengths or if the colloidal particles are covered by nonmagnetic or hardly 
polarizable spherical shells the dipolar energy at contact becomes comparable to ksT. In 
Subsec. |111 B| a previously developed density-functional theory has been applied to calculate 
the phase diagram including the temperature as a relevant thermodynamic variable. The 
overall effect of considering crystal structures different from fee on the position of the phase 
boundaries is rather small. The occurrence of a stable ferroelectric liquid phase within a 
certain parameter range is confirmed and in agreement with corresponding conclusions based 
on simulation data. Nonetheless it is conceivable that some trends in the behavior of the 
axis ratios, such as the density dependence of c/a in the hexc phase and the value of b/a in 
the bco phase, are not correctly described by DFT, which contains a number of uncontrolled 
approximations. However, at present no better theory for the quite demanding problem of 
freezing of dipolar fluids is available. 



APPENDIX: THE DIPOLAR LATTICE SUM 



The dipolar sum is only conditionally convergent, which means that its value depends on 
the order of summation, or, equivalently, on the sample shape. To clarify this often overlooked 
difficulty, here we explicitly perform the thermodynamic limit starting from finite lattices 
and letting the sample size diverge for fixed shape. We consider rotational ellipsoids with 
axis lengths kL along the z direction and L along the x and y directions. 
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1. Reduction to a single lattice sum 



The dipolar energy per particle is 



v dip 



^EE w ^ K - R ') = \ E h (^r) ^( R ^) (ai) 

R R'^R Ri27^0 ^ ' 

where the superscript L refers to a finite system and the lattice vectors R and R' run over 
the N sites inside the sample. In the second form one summation has been carried out 
so that R12 runs over a sample of doubled size, and the function h counts the number of 
occurrences of a given interparticle vector R 12 . If for large systems the discrete array of sites 
is approximated by a uniform distribution of the same density, the function h is given by 
the ratio of the intersection volume of two ellipsoids shifted by R12 relative to each other 
and the volume of one ellipsoid. The explicit result derived in Ref. |ll] is 



with 



3 / . - „ 1 



1/2 



/>i(#) = -| (sin 2 fl + — cos 2 fl ) (A3) 



and 



If 1 \ 3/2 

/i 3 (fl) = -(sin 2 £ + -cos 2 £j (A4) 

where 9 is the angle between R12 and the z axis. For rapidly decaying potentials in the 
thermodynamic limit L —>■ 00 h can be replaced by 1 because the difference h — 1 becomes 
appreciable only for large values of Ru comparable with L which, however, have a small 
weigth in Eq. ( |A1~1) due to the vanishing of w(R\2). But it is unclear whether this line of 
argument also holds for the slowly decaying dipolar potential. In order to check this we 
choose a cutoff radius R c beyond which one may approximate the summation in Eq. (|A1|) 
by an integral. In this limit the terms with i? 12 > R c become 

r 1 r 2L9(e) m 2 r 

2vrp / dcosO / drr 2 — P 2 (cos 6)h(0, -) 

J-l JRr ^ L 



'R, 

• 1 

,2 



2Tcpm / 0056^2 (cos 1 



1 



lng(9) + h 1 (e)g(e) + h 3 (e)g(d) 3 



+ 0(1). (A5) 



Here p is the density of dipoles, i^O 3 ?) — (3x 2 — l)/2 is the second Legendre polynomial, 
and the function g{6) = (sin 2 # + cos 2 O/k 2 ) -1 ^ 2 parametrizes the surface of the ellipsoid 



IIJ. Due to the specific forms of the functions hi and g the last two terms vanish so that 
indeed one obtains the same result if h = 1 is set from the beginning. The same is obviously 
true for the contributions from Ru < R c for L — > 00. These considerations justify that 
the average energy per particle may be replaced by the energy of the central particle. We 
emphasize that this argument hinges on the ellipsoidal shape of the sample. For example, 
by explicit calculations one can show that the replacement is not correct for a parallelepiped 
with general aspect ratio. 
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2. Shape dependence in the Ewald method 



The Ewald form of the dipolar lattice sum is often used without any discussion of the 
shape dependence of the original sum |41],|42|,|34[] . In the following we show that it actually 
corresponds to a specific choice of the sample shape and we derive the corrections that must 
be applied for other shapes. To this end the dipolar sum is constructed by superposition of 
two slightly shifted opposite point charge lattices. Hence we first recapitulate the derivation 
of the corresponding Ewald sum for the electrostatic potential 0(r) of a finite Bravais lat- 
tice of positive unit point charges plus an opposite uniform background charge. The Ewald 
method proceeds by rearranging the charge density p(r) into two contributions, employ- 
ing a suitable addition of zero: first pi(r) corresponding to a lattice of positive Gaussian 
charge distributions plus the negative background, second /?2( r ) corresponding to a lattice 
of negative Gaussians plus the positive point charges: 

p(r) = Pl (r) + p 2 (r) = $>< 0) (r - R) + pf (r - R)] (A6) 

R 



2 \ 3/2 ~ 

Pl 0) (r)=(-) e-^ 2 -le c (r) (A7) 



where 



and 

,(0 / • - * V 



7T 



tf>(T) = 5(r)- - e"— (A8) 



with c (r) = 1 if r is in the unit cell around the origin and c (r) = otherwise. The Fourier 
transform 



px(k) = / d 3 re- lkr Pl (r) (A9) 

of the first contribution is 

p L 1 (k)=pf\k)J2'e- lk - R . (A10) 



R 



The L dependence arises via the summation over the finite number of lattice vectors, indi- 
cated by the prime on the summation sign. The sum in Eq. ( |A10| ) is strongly peaked around 
the reciprocal lattice vectors G for large systems and approaches (2ti) 3 /V c ^ G 5(k — G) in 
the thermodynamic limit. The corresponding electrostatic potential is 



47T 



[ d 3 fce ikr pf(k)/fc 2 . (All) 



If the unit cell is chosen as the volume spanned by the basis vectors and as centered around 
r = one finds C (G) = for all G ^ and Pi(k = 0) = due to local charge neutrality. 
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For the second contribution p2(r) the potential is easily obtained by integration of the 
Poisson equation in spherical coordinates yielding 



R 



|r-R 



For L — > oo the total potential is 31 



, . . An 1 , r G 2 v-^ erfc(z/|r — R|) . ,. 

<^ r ) = V E ^ e exp( -^ } + E r _ R i + C'M (A13) 

c G^O R 

where C(z/) = — n/(i> 2 V^) is a constant stemming from the small k contributions in Eq. (|A11|) . 
In this form both sums are rapidly converging and 0(r) is independent of v. 

The potential of a corresponding arrangement of dipoles is now obtained from the su- 
perposition of two slightly shifted point charge lattices with opposite signs: 



,T f \ ,• m 



(r - d/2) - ^(r + d/2) + ^„,(r - d/2) - ^(r + d/2) 



d = dm. 

(A14) 



The potentials <p^ ni of the uniform background charges, which do not cancel each other in this 
limit, have been subtracted. Furthermore the hats indicate that the Coulomb potential 1/r 
must also be subtracted, because the field of the dipole at the origin must not be included. 
The energy of the central dipole follows as 



V. 



m d L | m 2 d 2 



di P 2 dz ^<%| r=0 2 dz 2 



r=0 



(A15) 



where the limit d — > generates the second derivative. The first term reproduces Eq. (0) for 
L — > oo (see Eq. ( |A13j )). The generalization to lattices with basis follows from summing the 
contributions from each sublattice. The second term is calculated using 

L ( (r=(0,0,r)) = p / d 3 r'-^— 



urn 



pi r2Lg(6) 

2tt P J dcos6 J dr'r' 2 {r 2 + r' 2 - 2rr' cos 6)~ 1/2 . (A16) 



By performing the inner integration and expanding in terms of r one finds 

9 xL ( r )L=n = - A7[ P dcos6P 2 {cos6)\ng{6)^ = -AixpD^k). 



dz 2 ^ umy y|r=0 r V3 



(A17) 



The expression in brackets is the depolarization factor D(k) of an ellipsoid with aspect 
ratio k All other terms are independent of k for L — > oo. Thus we have derived an 
explicit result for the shape dependence of the dipolar energy which is the same as for the 
homogeneous dipole density studied in Ref. [II]]. Since D{k — ► oo) = Eq. (|7]) is correct 
for a needle-shaped sample. In other sample shapes the dipolar energy can be lowered by 
domain formation. 



14 



REFERENCES 



[1] B. Groh and S. Dietrich, in New Approaches to Problems in Liquid State Theory, Vol. 529 
of NATO Science Series, edited by C. Caccamo, J.-P. Hansen, and G. Stell (Kluwer, 
Dordrecht, 1999), pp. 173-196, and references therein. 

[2] P. Teixeira, J. Tavares, and M. T. da Gama, J. Phys.: Condens. Matter 12, 411 (2000). 

[3] R. van Roij, Phys. Rev. Lett. 76, 3348 (1996). 

[4] J. Shelley, G. Patey, D. Levesque, and J. Weis, Phys. Rev. E 59, 3065 (1999). 
[5] D. Wei and G. Patey, Phys. Rev. Lett. 68, 2043 (1992). 

[6] J. Weis, D. Levesque, and G. Zarragoicoechea, Phys. Rev. Lett. 69, 913 (1992). 

[7] M. Stevens and G. Grest, Phys. Rev. E 51, 5976 (1995). 

[8] J. Weis and D. Levesque, Phys. Rev. Lett. 71, 2729 (1993). 

[9] M. Stevens and G. Grest, Phys. Rev. E 51, 5962 (1995). 
[10] B. Groh and S. Dietrich, Phys. Rev. Lett. 72, 2422 (1994). 
[11] B. Groh and S. Dietrich, Phys. Rev. E 50, 3814 (1994). 
[12] S. Klapp and F. Forstmann, J. Chem. Phys. 106, 9742 (1997). 
[13] R. Sear, Phys. Rev. Lett. 76, 2310 (1996). 

[14] M. Osipov, P. Teixeira, and M. T. da Gama, Phys. Rev. E 54, 2597 (1996). 

[15] B. Groh and S. Dietrich, Phys. Rev. E 54, 1687 (1996). 

[16] S. Klapp and F. Forstmann, Europhys. Lett. 38, 663 (1997). 

[17] S. Klapp and F. Forstmann, J. Chem. Phys. 109, 1062 (1998). 

[18] S. Klapp and G. Patey, J. Chem. Phys. 112, 10949 (2000). 

[19] R. Tao and J. Sun, Phys. Rev. Lett. 67, 398 (1991). 

[20] J. Martin, R. Anderson, and C. Tigges, J. Chem. Phys. 108, 3765 (1998). 

[21] T. Chen, R. Zitter, and R. Tao, Phys. Rev. Lett. 68, 2555 (1992). 

[22] U. Dassanayake, S. Fraden, and A. van Blaaderen, J. Chem. Phys. 112, 3851 (2000). 

[23] G. Gao and X. Zeng, Phys. Rev. E 61, 2188 (2000). 

[24] J. Weis and D. Levesque, Phys. Rev. E 48, 3728 (1993). 

[25] E. Lomba, J. Weis, and C. Tejero, Phys. Rev. E 58, 3426 (1998). 

[26] J. Fernandez and J. Alonso, Phys. Rev. B 62, 53 (2000). 

[27] B. Groh and S. Dietrich, Phys. Rev. E 53, 2509 (1996). 

[28] B. Groh and S. Dietrich, Phys. Rev. E 57, 4535 (1998). 

[29] A. Smith and N. Ashcroft, Phys. Rev. B 38, 12942 (1988). 

[30] B. Nijboer and F. de Wette, Physica 23, 309 (1957). 

[31] J. Slater, Insulators, Semiconductors, and Metals (McGraw-Hill, New York, 1967). 

[32] A. Denton and N. Ashcroft, Phys. Rev. A 39, 4701 (1989). 

[33] W. Curtin and N. Ashcroft, Phys. Rev. Lett. 56, 2775 (1986). 

[34] S. Gay, P. Beale, and J. Rainwater, J. Chem. Phys. 109, 6820 (1998). 

[35] R. Rosensweig, Ferrohydrodynamics (Cambridge University Press, Cambridge, 1985). 

[36] Electro -Rhcological Fluids, Magneto -Rheological Suspensions and Associated Technology, 

edited by W. A. Bullough (World Scientific, Singapore, 1996). 
[37] E. Dubois, R. Perzynski, F. Boue, and V. Cabuil, Langmuir 16, 5617 (2000). 
[38] E. Lemaire, G. Bossis, and Y. Grasselli, Langmuir 8, 2957 (1992). 
[39] J. Woestman and A. Widom, Phys. Rev. E 48, 1995 (1993). 

[40] L. Marshall, C. Zukoski, and J. Goodwin, J. Chem. Soc. Faraday Trans. I 85, 2785 
(1989). 



15 



[41] D. Adams and I. McDonald, Mol. Phys. 32, 931 (1976). 

[42] C. Kittel, Einfuhrung in die Festkorperphysik, 12th ed. (R. Oldenbourg, Munich, 1999). 



16 



FIGURES 

bco 




FIG. 1. Geometries of the first four lattice structures considered in the search for the energy 
minimum. For bco the situation at the fcc-bco transition is shown: when a small dipole moment is 
introduced in a nonpolar fee crystal, polarization along [110] is preferred and the crystal contracts 
in this direction. Thick and thin lines mark the conventional unit cells of fee and bco, respectively. 
Red particles lie on the face centers for fee and fco and at half height {z = c/2) for hexc. 
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hexab 




FIG. 2. The "hexab" crystal structure which arises from a distortion of a hexagonal close packed 
lattice polarized in a nearest-neighbor direction within the hexagonal plane. It is an orthorhombic 
structure with three additional basis particles shown in red (ti), green (T2), and blue (T3) with 
coordinates as given in the main text. Four orthorhombic unit cells are shown. The thin lines mark 
the hexagonal unit cell that is recovered for the special axis ratio c/a = l/\/3- 



18 



1 - 



C/) 

o 

CO 
■ — 

X 



0.9 



0.8 



0.7 



0.6 




m 

FIG. 3. Axis ratios of the bco phase of a Stockmayer solid at T = and p* = 1.24 as a function 
of the dipole moment. For m* = one has c/a = b/a = l/\/2 corresponding to fee. The diamonds 
denote simulation results of Gao and Zeng p| for p* = 1.24, m* = 2.5, and T* = 0.7. At m* ~ 4.25 
a continuous transition to bet takes place. Note that for m* < 1.95 the hexc structure is slightly 
more stable than bco, while for m* > 2.55 gas-solid phase separation occurs. Thus for this density 
bco is stable within the indicated window and metastable outside. 
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FIG. 4. Ground state phase diagram of the Stockmayer model. The solid (dotted) line denotes 
first (second) order transitions. The dashed line indicates the global minimum of the energy per 
particle over all densities. The unlabeled area to the left of this line is a two-phase region where 
an infinitely diluted gas coexists with the solid. The maximum of the pressure p = p 2 dU /dp in the 
depicted parameter range occurs in the lower right corner where p ~ 500e/cr 3 . Using parameter 
values for argon this corresponds to about 20 GPa which is accessible in a laboratory. 
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FIG. 5. Axis ratios of the bco structure of a Stockmayer solid at p* = 1.24 and T* = 0.7 
calculated from density-functional theory. Ferroelectric order sets in for m* > 0.67. The detailed 
behavior near this point could not be clarified due to numerical problems. 
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FIG. 6. Phase diagram of the Stockmayer fluid as calculated from density-functional theory for 
m* = 2. The ferroelectric solid has hexc structure. The dotted line indicates the gas-ferroelectric 
liquid-ferroelectric solid triple temperature. The dashed lines are the corresponding phase bound- 
aries if an fee solid is assumed. The unlabeled areas are two-phase regions. 
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FIG. 7. Ground state phase diagram of the dipolar soft sphere model for m* = 1. The meaning 
of the line styles is the same as in Fig. |j. The abbreviations for the various crystal structures are 
explained in the main text. The unlabeled areas are two-phase regions. The plots in the middle 
and at the bottom are magnifications of the regions around n = 12 and near the hard sphere limit 
n = oo, respectively. The phase diagram for other dipole moments m* can be inferred from the 
one given here by rescaling the reduced density according to p* /m* G ^ n ~ 3 \ In the limit p* — ► oo 
the upper bco phase turns into fee continuously. 
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FIG. 8. Interaction potential between parallel chains of dipolar soft spheres, with a longitudinal 
shift of half the intrachain particle distance a, for different values of the exponent n as a function 
of the chain separation r. The plot is for m* = 1; corresponding curves for other dipole moments 



can be obtained by rescaling the distance with m* 2 ^ n 3 ) and the energy with 



*2n/(n— 3) 
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